//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include <assert.h>
#include <math.h>

#include "ephi.hpp"

#define ABSURD_VALUE 9e9

Scene::Scene (ElectroDynamics& dynamics) : dynamics(dynamics)
{
  size_t i;
  for ( i = 0; i < 256; i++ )
      colors.push_back (Screen::color(255, i, 0));
  for ( i = 1; i < 256; i++ )
      colors.push_back (Screen::color(255 - i, 255, 0));
  for ( i = 1; i < 256; i++ )
      colors.push_back (Screen::color(0, 255 - i, 0));

#if 0
  Screen s(colors.size(), 30, 1.0, 1.0);
  for ( size_t i = 0; i < colors.size(); i++ )
    {
      for ( size_t j = 0; j < 30; j++ )
        {
          s.set_pixel (i, j, colors[i]);
        }
    }

  s.write ("grad.png");
#endif

  has_data = false;
  set_coloring ();
  step = 0.0;
  mbf = 0.0;
  bfadd = 0.0;
}

Scene::~Scene ()
{
}


Scene::Source::~Source ()
{
}

void Scene::calc (Screen& screen, ecalc_t ct, Source *src)
{
  vmap.init (screen.get_width(), screen.get_height());
  tmap.init (screen.get_width(), screen.get_height());
  smap.init (screen.get_width(), screen.get_height());

  max_bf = -ABSURD_VALUE;
  min_bf = ABSURD_VALUE;

  TaskManager tmgr;
  std::vector<calctask> tasks;
  tasks.resize (THREADS_CPUCOUNT);
  for ( size_t i = 0; i < tasks.size(); i++ )
    {
      tasks[i].part = (int)i;
      tasks[i].min_f = ABSURD_VALUE;
      tasks[i].max_f = -ABSURD_VALUE;
      tasks[i].scene = this;
      tasks[i].screen = &screen;
      tasks[i].ct = ct;
      tasks[i].src = src;
      tmgr.addTask (&tasks[i]);
    }

  while (1)
    {
      calctask *task = (calctask*)tmgr.run();
      if ( !task ) break;
      if ( task->min_f < min_bf )
          min_bf = task->min_f;
      if ( task->max_f > max_bf )
          max_bf = task->max_f;
    }

//  max_bf = 1.2; // XX
//  min_bf = 0.002;

  if ( min_bf < 0 )
    {
      bfadd = fabs(min_bf);
      Ephi::debug (Ephi::DEBUG_WARN,
          "Scene::calc: minimum field magnitude is negative (%g), adding %g to all field values\n",
          min_bf, bfadd);
    }
  else
    {
      bfadd = 0.0;
    }

  min_bf += bfadd;
  max_bf += bfadd;

  if ( max_bf / 100000 > min_bf ) // don't allow greater fields difference - in the interest of a good picture
      min_bf = max_bf / 100000;

  has_data = true;
  update_levels ();
}

//#define DDD
#ifdef DDD
#include <sys/time.h>
#endif

void Scene::run_calctask (calctask *task)
{
  Screen& screen = *task->screen;
  ecalc_t ct = task->ct;

  size_t sh = screen.get_height();

  for ( size_t y = 0; y < sh; y++ )
    {
      if ( (y % THREADS_CPUCOUNT) != task->part )
          continue;
#ifdef DDD
      struct timeval tv1, tv2;
      gettimeofday (&tv1, 0);
#endif
      for ( size_t x = 0; x < screen.get_width(); x++ )
        {
          vect3d p;
          screen.reverse (x, y, p);
          prec_t len = 0;
          vect3d f, tp;

          bool inContact = dynamics.inContact (p);
          switch (ct)
            {
              case CALC_BFIELD: {
                vect3d dummy;
                dynamics.getFields (p, f, dummy);
                vmap.set (x, y, f);
                len = f.length();
                screen.translate (screen.get_center() + f, tp.x, tp.y, tp.z);
                tp.z = 0;
                tp.normalize();
                tmap.set (x, y, tp);
                break;
              }
              case CALC_EFIELD: {
                vect3d dummy;
                dynamics.getFields (p, dummy, f);
                vmap.set (x, y, f);
                len = f.length();
                screen.translate (screen.get_center() + f, tp.x, tp.y, tp.z);
                tp.z = 0;
                tp.normalize();
                tmap.set (x, y, tp);
                break;
              }
              case CALC_POTENTIAL: {
                if ( !inContact )
                    len = dynamics.getPotential (p);
                break;
              }
              case CALC_SCALAR: {
                if ( !inContact )
                    len = task->src->getScalar(p);
                break;
              }
            }
    //          p.print("p");
    //          printf ("v=%g\n", (double)len);
          if ( inContact )
            {
              smap.set (x, y, ABSURD_VALUE);
            }
          else
            {
              smap.set (x, y, len);
              if ( task->max_f < len )
                  task->max_f = len;
              if ( task->min_f > len )
                  task->min_f = len;
            }
        }
#ifdef DDD
      gettimeofday (&tv2, 0);
      printf ("row %zu: %ld\n", y, tv2.tv_usec - tv1.tv_usec + 1000000 * (tv2.tv_sec - tv1.tv_sec));
#endif
    }
}

void Scene::set_coloring (size_t levels, bool use_log)
{
  if ( levels == 0 )
      levels = colors.size();
  this->levels = levels;
  this->use_log = use_log;
  update_levels ();
}

void Scene::update_levels ()
{
  if ( !has_data )
      return;
  if ( min_bf == 0 )
    {
      factor = 0;
      max_level = 0;
    }
  else
    {
      factor = prec_t_pow(max_bf / min_bf, 1.0/levels);
      max_level = prec_t_log(max_bf) / prec_t_log(factor);
    }
  range = max_bf - min_bf;
}

void Scene::render_map (Screen& screen)
{
  if ( !use_log )
      Ephi::debug (Ephi::DEBUG_INFO, "Scene: %u levels, min_sf=%.8g max_sf=%.8g, range=%.8g (%.2fx)\n",
          levels, prec2double(min_bf), prec2double(max_bf), prec2double(range), prec2double(max_bf/min_bf));
  else
      Ephi::debug (Ephi::DEBUG_INFO, "Scene: %u levels, min_sf=%.8g max_sf=%.8g, range=%.8g (%.2fx), logbase=%.8f\n",
          levels, prec2double(min_bf), prec2double(max_bf), prec2double(range), prec2double(max_bf/min_bf), prec2double(factor));

  size_t y;
  for ( y = 0; y < screen.get_height(); y++ )
    {
      for ( size_t x = 0; x < screen.get_width(); x++ )
        {
          prec_t m = smap[y * screen.get_width() + x];
          screen.set_pixel (x, y, cmap(m));
        }
    }
}

const Screen::color& Scene::cmap (prec_t bf)
{
  if ( range == 0 )
      return colors[colors.size() - 1];
  size_t r;
  bool is_absurd_value = bf == ABSURD_VALUE;
  bf += bfadd;
  if ( !is_absurd_value && bf > max_bf ) bf = max_bf;
  if ( is_absurd_value || bf < min_bf )
    {
      r = 0;
    }
  else if ( use_log )
    {
      if ( factor == 0 )
          return colors[colors.size() - 1];
      size_t level = (size_t)ceil(max_level - prec_t_log(bf)/prec_t_log(factor));
      if ( level == 0 ) // calculation errors at edges, I think, probably
          level = 1;
      else if ( level == levels + 1 )
          level = levels;
      level = levels - level + 1;
      assert(level <= levels);
      r = (size_t)prec2double(level * ((prec_t)(colors.size() - 1) / levels));
    }
  else
    {
      size_t level = (size_t)ceil(0.1 + ((bf - min_bf)/range) * levels);
      if ( level == levels + 1 )
          level = levels;
      assert(level <= levels);
      r = (size_t)prec2double(level * ((prec_t)(colors.size() - 1) / levels));
    }
  assert ( r < colors.size() );
  return colors[colors.size() - r - 1];
}

#define DS 0.0001

bool Scene::render_bline (Screen &screen, const vect3d& pos, const Screen::color& color, prec_t limit_dist,
        const vect3d& limit_pos, prec_t limit_pos_radius, bool backwards)
{
  prec_t factor = DS;
  prec_t fby6 = factor / 6.0;
  prec_t dist = 0;
  prec_t rps2 = limit_pos_radius * limit_pos_radius;
  prec_t ldist = 0;
  prec_t checkpoint_dist = screen.get_dx() / 3.0;
  vect3d P = pos;
  size_t old_x, old_y;
  screen.translate (P, old_x, old_y);
  enum { LEAVING, NEUTRAL, ARRIVING } status = LEAVING;
  prec_t last_mtoend = 0.0;
  prec_t sign = backwards ? -1 : 1;

  screen.set_pixel (P, color); // the initial pixel

  while ( true )
    {
      vect3d bf1, bf2, bf3, bf4;
      dynamics.getBField (P, bf1);
      bf1.normalize();
      dynamics.getBField (P + sign * 0.5 * factor * bf1, bf2);
      bf2.normalize();
      dynamics.getBField (P + sign * 0.5 * factor * bf2, bf3);
      bf3.normalize();
      dynamics.getBField (P + sign * 1.0 * factor * bf3, bf4);
      bf4.normalize();
      P = P + sign * fby6 * (bf1 + 2*bf2 + 2*bf3 + bf4);

      dist += factor;
      ldist += factor;

      if ( ldist > checkpoint_dist )
        {
          ldist = 0;
          if ( limit_dist != -1.0 && dist > limit_dist )
              break;
          switch (status)
            {
              case LEAVING: {
                if ( (limit_pos - P).magnitude() > rps2 )
                    status = NEUTRAL;
                break;
              }
              case NEUTRAL: {
                prec_t mtoend = (limit_pos - P).magnitude();
                if ( mtoend < rps2 )
                  {
                    status = ARRIVING;
                    last_mtoend = mtoend;
                  }
                break;
              }
              case ARRIVING: {
                prec_t mtoend = (limit_pos - P).magnitude();
                if ( mtoend > last_mtoend )
                    return true;
                last_mtoend = mtoend;
                break;
              }
            }
          size_t x, y;
          screen.translate (P, x, y);
          if ( x != old_x || y != old_y )
            {
              screen.set_pixel (P, color);
              old_x = x;
              old_y = y;
            }
        }
    }
  return false;
}

bool Scene::find_untouching_bline (vect3d& result, const vect3d& pos0, vect3d& pos1,
    prec_t limit_dist, prec_t limit_pos_radius)
{
  bool backwards = false;

  vect3d low = pos0;
  vect3d high = pos1;

  for ( size_t i = 0; i < 20; i++ )
    {
      if ( (high - low).magnitude() < DS * DS )
          break;
      vect3d pos = (high + low) * 0.5;
      vect3d limit_pos = pos;
      // pos.print("pos");

      prec_t factor = DS;
      prec_t fby6 = factor / 6.0;
      prec_t dist = 0;
      prec_t rps2 = limit_pos_radius * limit_pos_radius;
      prec_t ldist = 0;
      prec_t checkpoint_dist = DS * 3;
      vect3d P = pos;
      enum { LEAVING, NEUTRAL, ARRIVING } status = LEAVING;
      prec_t last_mtoend = 0.0;
      prec_t sign = backwards ? -1 : 1;

      bool is_touching = false;

      while ( true )
        {
          vect3d bf1, bf2, bf3, bf4;
          dynamics.getBField (P, bf1);
          bf1.normalize();
          dynamics.getBField (P + sign * 0.5 * factor * bf1, bf2);
          bf2.normalize();
          dynamics.getBField (P + sign * 0.5 * factor * bf2, bf3);
          bf3.normalize();
          dynamics.getBField (P + sign * 1.0 * factor * bf3, bf4);
          bf4.normalize();
          P = P + sign * fby6 * (bf1 + 2*bf2 + 2*bf3 + bf4);

          dist += factor;
          ldist += factor;

          if ( ldist > checkpoint_dist )
            {
              if ( limit_dist != -1.0 && dist > limit_dist )
                  break;
              bool do_break = false;
              switch (status)
                {
                  case LEAVING: {
                    if ( (limit_pos - P).magnitude() > rps2 )
                        status = NEUTRAL;
                    break;
                  }
                  case NEUTRAL: {
                    prec_t mtoend = (limit_pos - P).magnitude();
                    if ( mtoend < rps2 )
                      {
                        status = ARRIVING;
                        last_mtoend = mtoend;
                      }
                    break;
                  }
                  case ARRIVING: {
                    prec_t mtoend = (limit_pos - P).magnitude();
                    if ( mtoend > last_mtoend )
                        do_break = true;
                    last_mtoend = mtoend;
                    break;
                  }
                }
              if ( do_break ) break;

              if ( dynamics.inContact (P) )
                {
                  is_touching = true;
                  break;
                }
            }
        }

      if ( !is_touching )
        {
          high = pos;
        }
      else
        {
          low = pos;
        }
    }
  result = high;
  return true;
}

void Scene::calibrate_bline_step (const vect3d& pos, prec_t step)
{
  vect3d bf;
  dynamics.getBField (pos, bf);
  mbf = bf.length() * step;
  this->step = step;
}

void Scene::render_bline_range (Screen &screen, const vect3d& pos0, const vect3d& pos1,
    const Screen::color& color, prec_t first_factor, prec_t dist_limit,
    const vect3d& line_end, prec_t line_end_radius)
{
  vect3d normal = (pos1 - pos0).normal();
  vect3d pos = pos0;
  if ( step == 0.0 )
      calibrate_bline_step (pos0, (pos1 - pos0).length() / 10);
  prec_t acc = first_factor * mbf;
  prec_t sby10 = step / 10;
  prec_t ss = step / 100;
  vect3d dir = normal * ss;
  prec_t plane_d = -(normal * pos1);
  while ( pos * normal + plane_d < 0 )
    {
      if ( acc >= mbf )
        {
          char b[128];
          Ephi::debug (Ephi::DEBUG_INFO, "Rendering line at %s..\n", pos.sprint (b));
          if ( line_end_radius > 0 )
            {
              render_bline (screen, pos, Screen::BLUE, dist_limit, line_end, line_end_radius, false);
              render_bline (screen, pos, Screen::BLUE, dist_limit, line_end, line_end_radius, true);
            }
          else
            {
              if ( !render_bline (screen, pos, Screen::BLUE, dist_limit, pos, sby10, false) )
                  render_bline (screen, pos, Screen::BLUE, dist_limit, pos, sby10, true);
            }
          acc -= mbf;
        }

      pos = pos + dir;
      vect3d bf;
      dynamics.getBField (pos, bf);
      acc += bf.length() * ss;
    }
}

void Scene::render_lic (Screen &screen)
{
  if ( noise.init (screen.get_width(), screen.get_height()) )
    {
      // is this good for noise?!?
      prec_t *nend = noise + screen.get_width() * screen.get_height();
      for (prec_t *rit = noise; rit != nend; rit++ )
        {
          *rit = prec_t_drand();
          if ( *rit < 0.5 ) *rit = 0;
          else *rit = 1.0;
        }
    }

  TaskManager tmgr;
  std::vector<lictask> tasks;
  tasks.resize (THREADS_CPUCOUNT);
  for ( size_t i = 0; i < tasks.size(); i++ )
    {
      tasks[i].part = (int)i;
      tasks[i].scene = this;
      tasks[i].screen = &screen;
      tasks[i].noise = &noise;
      tmgr.addTask (&tasks[i]);
    }

  while ( tmgr.run() );
}

void Scene::run_lictask (lictask *task)
{
  Screen &screen = *task->screen;
  value_map<prec_t> &noise = *task->noise;

  size_t pc = screen.get_height() / THREADS_CPUCOUNT;
  size_t si = task->part * pc;
  size_t ei = si + pc;
  if ( task->part == THREADS_CPUCOUNT - 1 )
      ei = screen.get_height();

  for ( size_t y = si; y < ei; y++ )
    {
      for ( size_t x = 0; x < screen.get_width(); x++ )
        {
          vect3d P0((prec_t)x + 0.5, (prec_t)y + 0.5, 0);
          if ( smap.get(x, y) == ABSURD_VALUE )
            {
              screen.set_pixel (x, y, Screen::BLACK);
              continue;
            }
          prec_t total_h = 0;
          prec_t total_color = 0;
          lic_integrate (noise, P0, total_h, total_color, false);
          lic_integrate (noise, P0, total_h, total_color, true);
          prec_t f = total_color/total_h;
          Screen::color c = cmap(smap.get(x, y));
          if ( levels == 1 )
              c = Screen::WHITE;
          screen.set_pixel (x, y, Screen::color ((int)prec2double(f * c.r), (int)prec2double(f * c.g), (int)prec2double(f * c.b)));
        }
    }
}

#define HF_C 0.25
#define HF_D 0.75

prec_t Scene::hf(prec_t a, prec_t b, prec_t phi)
{
  prec_t ret = 0.25 * (
      b -
      a +
      (sin(b*HF_C) - sin (a*HF_C))/HF_C +
      (sin(b*HF_D + phi) - sin(a * HF_D + phi)) / HF_D +
      (sin(b*(HF_C-HF_D) - phi) - sin(a*(HF_C-HF_D)-phi)) / (2 * (HF_C - HF_D)) +
      (sin(b*(HF_C+HF_D) + phi) - sin(a*(HF_C+HF_D)+phi)) / (2 * (HF_C + HF_D))
    );
  return ret;
}

void Scene::render_lic_animation (Screen &screen, const char *name_prefix, int frames)
{
  value_map<prec_t> noise;
  noise.init (screen.get_width(), screen.get_height());

  // is this good for noise?!?
  prec_t *nend = noise + screen.get_width() * screen.get_height();
  for (prec_t *rit = noise; rit != nend; rit++ )
    {
      *rit = prec_t_drand();
      if ( *rit < 0.5 ) *rit = 0;
      else *rit = 1.0;
    }

  // this filter seems to help, together with external contrast stretching
  size_t bw = 4;
  size_t mx = bw * bw / 2;
  for ( size_t y = 0; y < noise.h - bw; y++ )
    {
      for ( size_t x = 0; x < noise.w - bw; x++ )
        {
          while (true)
            {
              size_t n = 0;
              for ( size_t i = 0; i < bw; i++ )
                {
                  for ( size_t j = 0; j < bw; j++ )
                    {
                      if ( noise.get(x + i, y + j) == 0.0 )
                          n++;
                    }
                }
              if ( n < mx ) break;
              noise.set (x + prec_t_lrand() % bw, y + prec_t_lrand() % bw, 1.0);
            }
        }
    }

  for ( int i = 0; i < frames; i++ )
    {
      prec_t phi = i * 2 * PREC_PI / frames;
      for ( size_t y = 0; y < screen.get_height(); y++ )
        {
          for ( size_t x = 0; x < screen.get_width(); x++ )
            {
              vect3d P0((prec_t)x + 0.5, (prec_t)y + 0.5, 0);
              if ( smap.get(x, y) == ABSURD_VALUE )
                {
                  screen.set_pixel (x, y, Screen::BLACK);
                  continue;
                }
              prec_t total_h = 0;
              prec_t total_color = 0;
              lic_integrate_hf (noise, P0, phi, total_h, total_color, false);
              lic_integrate_hf (noise, P0, phi, total_h, total_color, true);
              prec_t f = total_color/total_h;
              Screen::color c = cmap(smap.get(x, y));
              if ( levels == 1 )
                  c = Screen::WHITE;
              screen.set_pixel (x, y, Screen::color ((int)prec2double(f * c.r), (int)prec2double(f * c.g), (int)prec2double(f * c.b)));
            }
        }
      char buf[1024];
      sprintf (buf, "%s_%04u.bmp", name_prefix, i);
      screen.write (buf);
    }
}

#define LC 200
#define MAX_ADD_DS 10.0
#define CORRECTION 0.000001
#define PARALLEL 0.001

void Scene::lic_integrate (const value_map<prec_t>& noise, const vect3d& pos, prec_t& total_h,
    prec_t& total_c, bool neg)
{
  vect3d P = pos;
  prec_t add_ds = 0;
  prec_t add_h = 0;
  for ( size_t i = 0; i < LC; i++ )
    {
      prec_t ds = 42.0; // absurdly large value
      vect3d V = neg ? -tmap.get(P) : tmap.get(P);
      if ( fabs (V.y) > PARALLEL ) // V is not horizontal
        {
          prec_t bottom = (floor(P.y) - P.y) / V.y;
          if ( bottom >= 0 && bottom < ds )
              ds = bottom;

          prec_t top = (ceil(P.y) - P.y) / V.y;
          if ( top >= 0 && top < ds )
              ds = top;
        }
      if ( fabs (V.x) > PARALLEL ) // V is not vertical
        {
          prec_t left = (floor(P.x) - P.x) / V.x;
          if ( left >= 0 && left < ds )
              ds = left;

          prec_t right = (ceil(P.x) - P.x) / V.x;
          if ( right >= 0 && right < ds )
              ds = right;
        }
      ds += CORRECTION;

      prec_t hi = ds;
      add_ds += ds;
      add_h += hi;
      total_c += ds * noise.get(P);
      P = P + ds * V;
      if ( !noise.valid(P) )
          break;
      if ( add_ds > MAX_ADD_DS )
          break;
    }
  total_h += add_h;
}

void Scene::lic_integrate_hf (const value_map<prec_t>& noise, const vect3d& pos, prec_t phi, prec_t& total_h,
    prec_t& total_c, bool neg)
{
  vect3d P = pos;
  prec_t add_ds = 0;
  prec_t add_h = 0;
  for ( size_t i = 0; i < LC; i++ )
    {
      prec_t ds = 42.0; // absurdly large value
      vect3d V = neg ? -tmap.get(P) : tmap.get(P);
      if ( fabs (V.y) > PARALLEL ) // V is not horizontal
        {
          prec_t bottom = (floor(P.y) - P.y) / V.y;
          if ( bottom >= 0 && bottom < ds )
              ds = bottom;

          prec_t top = (ceil(P.y) - P.y) / V.y;
          if ( top >= 0 && top < ds )
              ds = top;
        }
      if ( fabs (V.x) > PARALLEL ) // V is not vertical
        {
          prec_t left = (floor(P.x) - P.x) / V.x;
          if ( left >= 0 && left < ds )
              ds = left;

          prec_t right = (ceil(P.x) - P.x) / V.x;
          if ( right >= 0 && right < ds )
              ds = right;
        }
      ds += CORRECTION;

      prec_t next_add_ds = add_ds + ds;
      prec_t hi;
      if ( neg )
          hi = hf (-next_add_ds, -add_ds, phi);
      else
          hi = hf (add_ds, next_add_ds, phi);
//      hi = ds;
      total_c += hi * noise.get(P);
      add_h += hi;
      add_ds = next_add_ds;
      P = P + ds * V;
      if ( !noise.valid(P) )
          break;
      if ( add_ds > MAX_ADD_DS )
          break;
    }
  total_h += add_h;
}

void Scene::renderGnuplotData (Screen& screen, const char *path)
{
  FILE *fp = fopen (path, "wb");
  assert (fp);
  float f;
  f = screen.get_height();
  fwrite (&f, sizeof(f), 1, fp);
  for ( size_t i = 0; i < screen.get_height(); i++ )
    {
      vect3d pos;
      screen.reverse (0, i, pos);
      prec_t x, y, z;
      screen.translate (pos, x, y, z);
      f = y;
      fwrite (&f, sizeof(f), 1, fp);
    }

  for ( size_t x = 0; x < screen.get_width(); x++ )
    {
      {
      vect3d pos;
      screen.reverse (x, 0, pos);
      prec_t x, y, z;
      screen.translate (pos, x, y, z);
      f = x;
      fwrite (&f, sizeof(f), 1, fp);
      }
      for ( size_t y = 0; y < screen.get_height(); y++ )
        {
          prec_t v = smap[y * screen.get_width() + x];
          f = v;
          if ( v == ABSURD_VALUE )
              f = max_bf;
          if ( f < 0 )
              f = 0;
          fwrite (&f, sizeof(f), 1, fp);
        }
    }
  fclose (fp);
}

